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Abstract 

This paper addresses the problem of consistent energy-based coupling of atomistic and con- 
tinuum models of materials, limited to zero-temperature statics of simple crystals. It has been 
widely recognized that the most practical coupled methods exhibit large errors on the atom- 
istic/continuum interface (which are often attributed to spurious forces called "ghost forces"). 
There are only few existing works that propose a coupling which is sufficiently accurate near the 
interface under certain limitations. In this paper a novel coupling that is free from "ghost forces" 
is proposed for a two-body interaction potential under the assumptions of either (i) one spatial 
dimension, or (ii) two spatial dimensions and piecewise affine finite elements for describing the 
continuum deformation. The performance of the proposed coupling is demonstrated with nu- 
merical experiments. The coupling strategy is based on judiciously defining the contributions 
of the atomistic bonds to the discrete and the continuum potential energy. The same method 
in one dimension has been independently developed and analyzed in [H. X. Li and M. Luskin, 
IMA J. Numer. Anal., to appear]. 

Keywords: atomistic model, consistent atomistic/continuum coupling, ghost force removal, atomistic bond 
contribution, multiscale method, finite element method 

AMS subject classification: 65N30, 70C20, 74G15, 74G65 

1 Introduction 

In many applications of solid mechanics, such as modeling cracks, structural defects, or nanoelec- 
tromechanical systems, the classical continuum description is not suitable, and it is required to 
utilize an atomistic description of materials. However, full atomistic simulations are prohibitively 
expensive; hence there is a need for efficient numerical methods that couple a continuum descrip- 
tion of a material in the region where material deformation is smooth and an atomistic description 
where the variations of the deformation gradient are large. 

All atomistic/continuum (A/C) coupling methods can be divided into two categories: the 
energy-based coupling and the force-based coupling. Energy-based coupling [HI El E21 ED ESI 
IM[ [35] consists of composing the energy of the system depending on both discrete and continuum 
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deformations, and driving the system to a stable equilibrium. A major challenge for energy-based 
A/C coupling is the presence of ghost forces — the spurious forces that create an error in the solution 
near the A/C interface that cannot be reduced by enlarging the atomistic region or refining the 
computational mesh in the continuum region. 

In this paper the term consistent will be used as a synonym for the absence of ghost forces. It 
has been recently proved in two dimensions (2D) under some rather general technical assumptions, 
that absence of ghost forces implies first-order consistency |22j . In addition, in [23] it is shown 
directly that the proposed method exhibits a first-order convergence rate. 

In one dimension (ID) the problem of consistent zero-temperature static energy-based coupling 
is easy to solve [HJ \T7\ 12 lj . However, in higher dimensions this is a true challenge. In the case 
of nearest neighbor interaction there are methods that are free from ghost forces (the coupling of 
length scales method [3] is an example of such a method), but in a general situation all existing 
coupling methods introduce a certain interfacial error. 

Several methods have been proposed to address the problem of consistent energy-based coupling, 
including the quasinonlocal (QNL) quasicontinuum method [31J and the geometrically consistent 
scheme (GCS) [8j [21], which may be regarded as a generalization of QNL. The latter method 
consists of two steps: first, passing from general finite range interaction to a nearest neighbor 
interaction in a padding region between the atomistic and continuum region, and second, passing 
from the nearest neighbor atomistic model to the continuum model. E, Lu, and Yang showed that 
the latter method exhibits no ghost forces in ID and no ghost forces in higher dimensions in the 
case of a straight (planar) interface with no corners [8j. 

Another noteworthy effort to minimize the ghost forces is a work of Klein and Zimmerman 
|12j who considered a problem of coupling in arbitrarily overlapping atomistic and continuum 
regions. Assuming a two-body interatomic potential, they proposed a method of computing the 
contributions of the atomistic bonds in the overlapping region by numerically minimizing the ghost 
forces using the least squares technique. 

The analysis of energy-based methods reveals the following advantages of consistent energy- 
based methods in ID: (i) their error can be efficiently controlled as opposed to a finite error of 
nonconsistent methods [U \21\ [23] , and (ii) their region of stable deformations essentially coincides 
with the one of the atomistic model, whereas inconsistent methods significantly underpredict the 
critical strain [H [7] . 

In the force-based methods [TOj, [131 EB1 1211 EQ] , the equilibrium is achieved by computing 
the generalized forces for each degree of freedom (associated with the atoms and the finite element 
nodes) and driving them to zero. The force-based methods are the most commonly used A/C 
coupling methods, as most of them exhibit no ghost forces by construction. However, there are 
difficulties arising from the force field being nonconservative [51 El [20] . 

In summary, consistent energy-based methods have good accuracy and stability properties; 
however, the existing energy-based methods are consistent only in a few special cases. The force- 
based coupling is the most popular alternative to the energy-based coupling, but the force-based 
methods are not sufficiently well understood at present. 

In the present paper a new consistent energy-based A/C coupling is proposed. The coupling 
is based on judiciously treating the atomic bonds near the A/C interface by consistently defining 
their contributions to the discrete and the continuum energy. Two variants of the method are 
formulated; the first one couples the atomistic and the continuum energies through the atoms 
inside the continuum region, and the second one performs coupling only through the interface. 
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The scope of the present paper is limited to a two-body interaction potential, one or two spatial 
dimensions, and piecewise linear finite element discretization of the continuum deformation (the 
latter assumption is only required for the two-dimensional case). No restrictions on the finite 
element mesh (except that its nodes are positioned at the lattice sites) are made. 

We note that essentially the same method in one-dimensional setting has been independently 
proposed and analyzed by Li and Luskin [TC] . 

The paper is organized as follows: Section [2] introduces the problem of A/C coupling and 
briefly discusses major difficulties. In section [3] the proposed methods are presented in detail in 
a one-dimensional setting and then extended to a two-dimensional case in sections [4] and [5j The 
numerical experiments illustrating performance of the proposed methods are presented in section 
[6j The paper ends with a discussion of the results (section [7]) and concluding remarks (section [8]) . 



2 A/C coupling in ID 

In order to give a more vivid illustration of the proposed methods, we will start with a one- 
dimensional problem formulation followed by a detailed presentation of the proposed methods in 
ID (section [3]) . We will then show how to extend the proposed methods to 2D in sections [|and[5| 



2.1 Atomistic model 

Consider a one-dimensional atomistic material described by positions of atoms in the reference 
configuration as Xi = i (i G Z). For simplicity we set the lattice parameter e = 1. Let I C Z be a 
set of indices of atoms present in the atomistic material. Atoms may displace from their reference 
positions X{ to positions yi, and their displacements are Ui = yi — Xj. The generic linear space 
containing Xi, yt, and Ui will be denoted as IA = {g : X — > M}. For elements of U introduce a finite 
difference operator D r as follows: 

{D r g)i = D r gi = g i+r - g r (for r G Z). 

Note that if i + r ^ X, then D r gi is undefined; therefore, strictly speaking, D r should not be 
considered as an operator U U. 

The interaction of the atoms is described by the interaction potential ip and the cut-off radius 
R, yielding the total interaction energy 

E{ y )= v{y 3 -yi)= <P{D r yi). (2.1) 

l<j-i<R l<r<R 

Here and in what follows we assume that yj > yi for j > i. In practice, the proximity of atoms is 
measured in the physical domain, i.e., the atoms i and j contribute a nonzero interaction energy if 
\yj — Vi\ < R- However, for the purpose of studying consistency of the method, this treatment is 
essentially the same as \j — i\ < R, since these criteria are equivalent on the uniform deformation 
Vi = Fxi. 

We consider the Dirichlet-type boundary conditions; namely, we fix the positions yi of atoms 
j G Id near the boundary: 

V i = Fxi (i G 2d), (2.2) 
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Figure 2.1: An illustration of the assumption (2.3): Any unconstrained atom % G has a full 



set of R left and R right neighbors inside I. The interaction of the leftmost unconstrained atom 
with the three left and the three right neighbors is displayed, constrained atoms i G Id are shown 
as smaller points, and the illustration is for R = 3. 

where F is an arbitrary deformation tensor. To avoid boundary effects, we require that every 
unconstrained atom i has a full set of neighbors: 

(Vi G X \ Z D ) (Vr G Z : |r| < -R) i + r G X. (2.3) 



(See Figure 2.1 for illustration.) 

The boundary conditions span the manifold of admissible deformations and the linear space of 
test functions 

= {|/£W: Vi = Fxi for all i G Xd}, Uq = {u G U : Uj = for all i G Id}. 
Under the assumption ( |2.3[ ) , the elements of satisfy 

^ D r vi = for all u G U , r = 1, . . . (2.4) 

Compute the variation of E(y): 

E'{y;v)= ^2 V 7 '{D r yi)D r Vi. 

l<r<R 

The equilibrium equations of the atomistic material under the external force fi in variational form 
can be written as 

find y G U D : E'{y; v) = for all u G W . (2.5) 

These equations admit the solution described by a uniform strain y = Fx, as is shown by the 
following computation: 

E'(Fx;v)= tp'(rF)D rVi = ^ <p'(rF)( ^ D r vA = for all veU , (2.6) 

Kr<R 



which is due to (2.4) 



The problem (2.5), although discrete, is usually too large to handle on a computer. Therefore, 
its approximations with reduced degrees of freedom are used. Normally, to have a converging 
numerical method with the energy E(y), one must have E'(Fx; v ) = for all v G Uq. This relation 
is sometimes called the "patch test", a term borrowed from the theory of finite elements to describe 
a necessary condition for nonconforming elements to converge [32]. We call such approximations 
E{y) consistent. It should be noted that it has been recently proved in 2D, under rather general 
assumptions on the mesh, the atomistic region, and the coupling mechanism, that the absence of 
ghost forces implies first-order consistency |22j . 
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2.2 Continuum model 



If the deformation gradient y 



t+i 



?/j is smooth in a neighborhood of some domain fi c , the exact 



atomistic energy (2.1) can be approximated with the continuum one 

E c {y) 



R 

£ 

r=l 



(/?(V r y)dx, 



(2.7) 



where V r = r-^ and y € VF 1,00 (f2 c ) is a continuum approximation to the discrete deformation (i.e. 



Vi ~ y(i))- The formula (2.7) is essentially the Cauchy-Born rule [2], where we approximate the 
interaction energy using the energy density Ylr=i ^C^rllix))- 
Compute the variation of E c (y): 



E' c (Fx,v) = Y^ [ f'(rF)V r vdx= (^2(p'(rF)r\ f 

r=l J ^ r=l ' X 



dx 



v{x)dx. 



If the deformation at the boundary of the continuum region is fixed, i.e., the admissible deformations 
v satisfy v\gn c = 0, then, as follows from application of the Green's formula to the above calculation, 
the uniform strain y(x) = Fx is also an equilibrium in the continuum model. 

2.3 Problem formulation 

Consider an atomistic material which in its reference configuration occupies the region f2 = (-N — 
R, N + R). The material will be treated continuously in the continuum region £l c = (0, N) and 
discretely in the atomistic region Q a = Q \ £l c (here • denotes the closure of a set). The A/C 
interface is T := dQ c = {0,iV}. Define the atomistic lattice I = (idZ and the atoms within 
the atomistic region X a = fi a n 7L. The set of atoms involved in formulation of the Dirichlet-type 
boundary conditions (2.2) is chosen as Zd = {i £ Z : N < \i\ < N + R}. For these I and Id, the 



relation (2.3) holds, and hence (2.4) also holds. 



Note that from an algorithmic point of view, fixing the positions of atoms i > N + 1 (i.e., near 
the right boundary) is redundant, since for the continuum model in the region (0, N) it is sufficient 
to fix the position only for the atom i = N. Nevertheless, to compare the approximate model with 
the fully atomistic one, it is convenient to keep these atoms fixed, as we will see below. 

It should also be noted that the regions were chosen in this way only for ease of visualization, and 
the discussion below is valid for a more general family of regions £l a , fi c , and Ijy. After presenting 
and discussing the proposed methods in ID, we will extend the method to the two-dimensional case 
where we will assume a general form of the regions. 

The deformation of the material will thus be defined on X, U T U fL: 



U = {y : X a U fl c 
Uo = {u £ U : Ui 



for all i £l D }. 



A typical element of IA is illustrated in Figure 2.2 Here and in what follows we interchangeably 
use the notation y% = y{i) = y{x{) for the point values of y G U. 
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Figure 2.2: A typical element y oiU, for N = 10 and R = 3. The smaller points correspond to 2d- 
2.4 Quasicontinuum method 

In this subsection we formulate the energy-based quasicontinuum (QCE) method [33J, show that it 
has ghost forces, and describe the existing energy-based strategies of removing them. The presence 
of ghost forces has been widely discussed in the literature [U EJ El H3 EI]; hence the present 
discussion will be very brief and will only serve to contrast the ideas of the proposed method with 
the existing methods. 

QCE can be defined as follows. Write the atomistic energy in the form 



l<[r[<R 

where we define <j){z) = ip(\z\), and interpret the expression in the parentheses as the energy 
associated with a particular atom i £ I. Then substitute the atomistic energy associated with 



atoms i ^ f2 a U T in the expression (2.8) by the continuum energy: 



E &)= E i\ E myi)) + Y,{\ E tfrvi)) 

ieXaUr ^ i+reX ' i&X c ^ i+r&X ' 

l<\r\<R l<k|<i? 

« E (\ E w^ + eQ/? E t(y r y)dx) 

ieXaUr ^ i+r&X ' ieX c ^ 1 2 i+reX ' 

l<|r|<R l<\r\<R 

E (\ E HD ryi ))+ir 2 E ^(V,y)dx=:^ ce (y), 



i6XaUr v i+reX 7 2 i+re2 

l<\r\<R l<\r\<R 



where the interval (i — ^, i + ^) is the "effective volume" of an atom i £ I c . 

To show that such a coupling has ghost forces, we let R = 2 and show that the variation of 
the energy (E qce )'(y;v) for a uniform deformation y = Fx does not vanish. Indeed, omitting the 
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details of straightforward but tedious calculations, one gets 

(E^)'(Fx;v)= E 4>'{rF)D r v^\ + l - £ <j>' (rF) F * V r vdx 

ieZ a ^ i+rel ' i+rel 3 

r=±l,±2 r=±l,±2 

//EA«l+«0, // OE ,N«2+Vl+«t)+«-l 

= ¥>(^) — 2 — + ( ^( 2F ) 2 

+ <//(F) + <//(2F) 2v{\) + (boundary terms), 

where the boundary terms are concentrated near i = N and are not relevant to the A/C coupling. 
This can be further simplified, assuming a piecewise affine interpolation in f2 c (then v (n) = ^4^ ): 

(^ cc )'(Fx;?;) = ip'(2F) V2 - Vl ~ - + + (boundary terms), 

which is clearly nonzero. 

The interpretation of presence of the nonzero (ghost) force, for instance, on atom i = — 1 can be 
as follows. The interaction of atoms i = — 1 and i = 1 is computed two times in a different manner: 
once according to the continuum strain 2y'(x\) and once according to the exact strain y\ — y~\. 
The first computation contributes the force only to the continuum region, which causes a loss of 
balance of forces on the atom i = — 1. 

The energy-based approach to removing the ghost forces is to modify the energy E qce (y) in such 
a way that the uniform deformation y(x) = Fx is the solution to the equations. In the context of the 
quasicontinuum method, this was first done in [31] for the second nearest neighbor (i? = 2) and then 
generalized in [8] to longer interactions. However, when generalizing these methods from dimension 
one to higher dimensions, additional difficulties arise, for instance, related to transition between so- 
called element-based and atom-based summation rules near interface corners (i.e., interface edges 
and vertices in three dimensions) [8]. 



3 Consistent A/C coupling in ID 

We propose a strategy of constructing a consistent coupling of the atomistic and the continuum 
models. The strategy is based on splitting the energy of atomic bonds into atomistic and continuum 
contributions. 



3.1 Exact and continuum contributions of a bond 

We start with introducing some preliminary terms and definitions. The term bond between atoms 
i £ Z and % + r G Z will refer to an open interval b = (i,i + r). Introduce the set of all bonds in the 
atomistic system 

B := {(i, i + r) : 1 < r < R, i G X, i + r G X}. 
We denote the potential energy of a bond b = (i, i + r) as 



(3.1) 
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and call it the exact contribution of the bond (i, i + r) to the potential energy. Thus, the exact 
atomistic potential energy (2.1) can be written as 

E{y) = Y,^{y). (3.2) 

beB 

For a bond b = (i,i + r) fully contained in Q c we define its contribution to the continuum energy 
or, in short, continuum contribution, as 



c b (y) := ~ J (p(V r y)dx (if b C fi c ). 



Note that we later extend the definition of c b (y) to all bonds in (3.5). 

The following proposition states that the variation of the exact contribution of a bond b coincides 
with the variation of the continuum contribution of 6 on a uniform deformation y = Fx. 

Proposition 3.1. 

c' b (Fx; v) = e' b (Fx; v) for any F > 0, b G B, v £ U. 
Proof. The validity of this proposition is verified by the following straightforward calculation: 

c[^ i+r) (Fx;v) = \ j <p'(rF)r-^dx = <p'{rF)v\ x £ r = Lp'(rF){v i+r - Vi), 

e[ iji+r ^(Fx;v) = Lp'(rF)D r Vi = ^{rF){v i+r - Vi). 

□ 



3.2 Method of combining exact and continuum contributions 



We now define the proposed A/C coupling method motivated by Proposition 3.1 

E ecc (y) := + E °b(y) =■ E ? c (y) + £ c ccc (y), 



(3.3) 



beB 

b<£ n c 



beB 
bcttc 



which is obtained by substituting e b (y) with c b (y) in (3.2) for the bonds b fully contained in 
the continuum region. We hence name it the method of combining the exact and the continuum 
contributions (ECC) of the bonds (hereinafter referred to as the ECC method). 



Proposition 3.2. The ECC method (3.3) is consistent, i.e., 

(E ecc )'(Fx;v) = 0. 



Proof. The consistency of ( |3.3[ ) follows directly from Proposition 3.1 and the identity (2.6): 
(E ecc )'(Fx;v) = e b (Fx;v) + £ c' b (Fx;v) = ^ e' b (Fx; v) = E'(Fx;v) = 0. 



beB 
b<£n c 



beB 

bCSl c 



beB 



□ 
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The remainder of this subsection will be devoted to showing that E ecc (y) can be computed 
efficiently, i.e., without a need to go through all the bonds b G f2 c - More precisely, we will show the 
validity of the following proposition. 

Proposition 3.3. The energy of the ECC method can be written as 

E ecc (y) = ( e ^) - c b(y» + E M, (3-4) 

beB 
b<£Q c 

where the definition of cj,(y) is extended to all bonds b as 

c b{y) '■=- I <p{V r y)dx. (3.5) 
r J 

f2 c nfc 

According to this proposition, assembling the energy can be done as follows: We go through all 
the atoms in the atomistic region, computing the energy of their interaction with all other atoms 
and subtracting continuum contribution of the energy (3.5) for the bonds crossing the interface 
(i.e., for the bonds b Q c such that ft c D b is nonempty). The energy thus computed is then added 
to the continuum energy E c (y). Seen in this way, the method can be efficiently implemented with 
a single loop over the atoms within the atomistic region only. 

Before we commence with a proof of Proposition 3.3, let us formulate, without a proof, the 
following result which is trivial in ID but will be crucial in extending the method to 2D. 

Lemma 3.4 (one-dimensional bond-density lemma). Almost each point x 6R is covered by exactly 
r (r G Z + ) bonds of the form (i, i + r) (i £ Z), i.e., 

E-X(i,i+r)(x) = 1 for all r G Z + , (3.6) 
r a.e. 

where x»( x ) * s ^ e characteristic function of a set. 

Proof of Proposition 3.3. First, fix r G Z + and notice that, as long as (2.3) holds, any bond b = 
(i, i + r) having a nonzero intersection with fi c belongs to B. Hence we can replace the sum over 
all i € Z in (3.6) by the sum over (i, i + r) G B: 

EX(ii+r)( x ) = 1 for all x G fi c , r G Z + . (3.7) 
v ' ' a.e. 

(i,i+r)eB 

Hence the continuum contribution of all the bonds in the system is 

^Z c b(v)=^Z\ I <f{V r y)dx = ~ / Xb(x)<f(V r y)dx 

tp(V r y)dx = E c (y). 
n 

Therefore we can write E^ cc (y) (cf. ( |3.3[ )) in the form 

£c CC (y)= c b (y) = E c (y) - ^c b (y), 
beB beB 

bcQ c b<£fl c 



and hence (3.4) follows. □ 
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Remark 3.1. The same coupling has been independently proposed and analyzed by Li and Luskin 
in a one- dimensional discrete setting. In I16\j they introduce an extension of the QNL method ]3T$ 
to the finite-range potential (i.e., arbitrary R), which is essentially equivalent to the ECC method 
of the present work, and they analyze it in a linearized discrete setting without coarsening (i.e., 
for functions y(x) which are piecewise affine on each interval i — 1 < x < i, i = 1, 2, . . . , N). 
In particular, they prove that (i) under certain assumptions the uniform deformation is a stable 
equilibrium, and (ii) the coupled method converges to the exact atomistic solution. 



3.3 Method of combining atomistic and continuum contributions 

The method proposed in the previous subsection consists of treating each bond that does not fully 
lie in the continuum region atomistically and modifying the continuum energy by subtracting the 
corresponding contribution from the continuum energy. Since the positions of atoms which are 
strictly inside Q c enter the expression of the total energy, the coupling between regions becomes 
nonlocal. For implementation it could be more preferable to have a coupling through the interface 
only. In this subsection we present such a method. The details of construction of this method, 
however, are somewhat involved and will be needed only to better understand a similar method in 
2D (section^. 

The calculation (3.8) indicates that the method with the exact continuum energy E c {y) should 
contain all the continuum contributions c&(y). Then we should derive an atomistic contribution 
a (i,i+r) (y) which would consistently balance c^_j_ r ) 

To do that, define the operator D^y for uj C O (uj ^ 0) in the following way. If uj = U m =i(^ m ' rm ) 
is a union of nonintersecting intervals, then 

M 

Du>y ■= ^2 (y(r m ) - y(l m )), 

m=l 

where \uj\ is the total length of uj. D u y has the following properties: 

1. rKD^y approximates the derivative of y, in particular, ^D w (Fx) = F. 

2. D( iji+r) y = D r yi. 

3. If uj = oj\ U u)2 and ui\ n u)2 = 0, then D^y = D Ul y + D U2 y. 

4. J v'(x)dx = D u v. 

Using this operator define the atomistic contribution of a bond b = (i,i + r) as 

a b (y) := (3.9) 



Compute the variations of the atomistic and the continuum (cf. ( |3.5[ )) contributions of a bond 
+ r) = b E B, using properties 4 and 1 of D^: 



1 f* (^Lv 

c' b (Fx;v) = - / (p'(rF)r— dx =tp'(rF) D bn n v and 
r J ax 



bnn c 



a' b (Fx; v) = (V (rF) — ^— D W ) =^'{rF) D bnn& v, (3.10) 
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and therefore, using properties 3 and 2 of D w , one can see that 

c' b (Fx; v) + a' b (Fx; v) = <p'(rF)D b v = <p'(rF)D r Vi = e' b (Fx; v) 
This immediately implies that if we define the energy 

E^{y) :=Y^a b (y)+Y. c ^y) = : 

beB beB 

then the coupling given by E licc (y) is consistent: 

(E acc )'(Fx;v) = ^e' b {Fx;v) = E'(Fx;v) = 0. 

beB 

We will call it the atomistic and continuum contributions (ACC) method. 
Due to (3.8) the method can be written as 

E^ c (y) = Y J ^b(y)+E c (y)- 

beB 



(3.11) 



Note that in (3.11 ) the sum over b G B can be effectively changed to the sum over (b E £>, b <£_ O c ), 
since if b C O c , then a b (y) = 0. 

The prominent feature of the method is that the atomistic part of the energy depends only on 
the deformation in the atomistic region and on the interface. This may be more convenient for 
implementation (see section 5.2) and can potentially help in parallelization of the method. 

We conclude this discussion by presenting another version of the atomistic contribution a b (y). 
if b n n c is a union of nonintersecting intervals (( m , f]rri) (1 < tti ^ then we define 



M 

~ I \ \ " Vm — ?m / 

oft(y) : = 2^ — : — H 



m=l 



<p[r- 

Vm ?m 



Its variation 



M 



a b {Fx;v) = ^ - 



Vm £,m v Vm v £n 



M 



'-ip (rF)r- 



m=l 



lm sra 



v'irF)^ 



Vm 



m=l 



coincides with the variation a' b (Fx;v) (cf. (3.10)), and hence a b (y) can be substituted by a b {y) in 
the approximation (3.11): 

E* cc (y) = ^2~a b (y)+E c (y). 
beB 

The approximation E &cc {y) has a potentially smaller number of coupled terms when the bonds 
cross the interface several times (this is more likely to happen in many dimensions), which may be 
easier for implementation. In the present work the ACC method was implemented for a convex 
atomistic region, in which case both methods coincide. 



CONSISTENT ATOMISTIC/CONTINUUM COUPLING 



12 



4 Extension of ECC to 2D 

Coupling atomistic and continuum models of materials becomes harder in more than ID. For 
instance, the quasinonlocal quasicontinuum method |3 1 j and its generalizations [SJ [21] suffer from 
ghost forces near the corners of an A/C interface. 

In this section we will extend the ECC method to 2D. The method will be consistent by 
construction, with no restrictions on the mesh (except that its nodes coincide with lattice sites). 

In the two-dimensional case the reference configuration is usually described by a uniform lattice 
Xi = f\i (where i S Z 2 , X{ E M 2 , A £ M 2x2 ), and the deformed configuration yi is often considered 
as given by some mapping Y : M 2 — > M?. However, in this paper we adopt a slightly different 
point of view where we set Xi = i and the uniform deformation tensor A will be accounted for as 
Hi = Y(Axi) (this is done, for instance, in [E]). This point of view is actually closer to computer 
implementation, where the lattice is indexed with integers rather than real numbers. We stress 
that this is not a limitation of the proposed methods; in fact the numerical examples (section [6]) 

will be presented for a hexagonal lattice with A = ^ ' 
4.1 Preliminaries 

Consider an open bounded region Q, C M 2 , the atomistic and the continuum regions f2 a C Q, 
J7 C = Q \ il a (both open), and the A/C interface T = <9f2 c . We also consider the atomistic lattice 
X = Q n Z 2 , the atoms within the atomistic region X a = f2 a n Z 2 , and the set of atoms involved 
in posing the Dirichlet-type boundary conditions Xq Clfl f2 a . To avoid boundary effects, we will 



make additional assumptions on X and Xd after we introduce atomistic interaction (in section 4.2). 

We further assume that J7 C is a polygon with vertices coinciding with some lattice sites Xi. We 
will consider the fully discrete case; i.e., we introduce a triangulation T of £l c with triangles T £ T 
whose vertices are also positioned at the lattice sites x%. The space of continuous piecewise affine 
finite elements on O c is denoted as Vi(J~). This setting is exactly the same as the one of the QCE 



method [33J. However, the difference between the QCE and the proposed methods will be in the 



way the atomistic and continuum energies are coupled (cf. section 4.5 for more details on relation 
between ECC and QCE). 

To formulate the proposed coupling in 2D, we introduce the following supplementary notations. 
For u : X — >• K, i, i + r G X, define the discrete differentiation 

(D r y)i = D r yi = y i+r - yi. 

Define a bond {x\,X2) between two points x\,xi 6 Z 2 as an interval 

(a*, x 2 ) := {(1 - A)a?i + Xx 2 : A € (0, 1)}. 

Define the averaging over a bond b = (x\, x 2 ) of a piecewise continuous function / : b — > M as 



1 

/ /(x)db := f /((l - A)*! + 



\x 2 )d\. 



The following property then holds: 

/ V r /db = A./i for any / G V\(T), 
J (i,i+r) 
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where the directional derivative is defined as 

Vrf(x) := lim(/(x + er) - f{x)). (4.1) 

If / is smooth, then V r / = r • V/. It should be noted that V/ of a function / G V\{T) may 
be undefined on a bond which fully or partly lies on an edge of some T G T ■ Nevertheless, the 
directional derivative V r / is piecewise continuous (more precisely, piecewise constant) on any bond 
b = (i, i + r) C f2 c with the same direction vector r. 

For vector-valued functions v the bond averages and the directional derivatives are defined 
componentwise in the same manner. 

4.2 Atomistic model and continuum approximation 

We assume that atomistic interaction is given by a set of neighbors TZ C Z 2 \ {0} and a two-body 
potential ip. For simplicity of notations we denote 

(j){z) := <p(\z\) for z G I 2 . (4.2) 

Denote the collection of all bonds in the system by 

B : = {(*, i + r): r G K, iel, i + r G 1} 

and the exact contribution of a bond b = (i,i + r) G B under the deformation y as 

e b(y) ■= (p(DrUi)- (4.3) 

It should be understood that in this formula r and i essentially depend on b, but to simplify the 
notations we will avoid writing r& or %. 

The energy of the atomistic model then reads 

E{y) = Y,Zb{y), (4.4) 
beB 

and its continuum approximation in Vt c based on the Cauchy-Born rule is 

E dy)= [Y,^y)dn (4.5) 

= E E l T l <M v ^It) (for y G Vi(T)). (4.6) 



TeT r&l 



To avoid boundary effects entering our analysis, similarly to the assumption (2.3), we assume 



(W G I \ 2d) (Vr G ft) i + rel and (4.7) 
(VrGft) (ViGZ 2 : (y + r)nO^/l) M + r G Z. (4.8) 

Practically, these assumptions mean that we have enough atoms Id whose position we fix so 
that each free atom near the boundary of (and hence each free atom in f2 c C f2) has enough 
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neighbors to interact with. It can be achieved by choosing a region near the boundary of £1 
such that H f2 c = and 

dist(J7 \ S1d) dQ) > max \r\ 



(see illustration on Figure 4.1 ). Then 2d = f^D H Z will satisfy (4.7) and (4.8). 

The spaces of deformations and displacements are then defined similarly to the one-dimensional 
case: 

U = {y:l eL uU~ c ^R 2 : y\ Uc G Pi(T)}, 
Uq = {u £U : Ui = for all i G 2d)}. 



With the assumptions (4.7) and ( |4.8| ), one can now show that the uniform deformation y = Fx is 
an equilibrium of the atomistic energy (4.4). 

Proposition 4.1. 



E'(Fx\ v) = for all F e M 2x2 , u G 



Proof. Compute 



e' (w+r) (Fx;t;) = </>'(Fr) • D r v u 
where, due to (4.2), 4>'{z) = y'd^Drfr f° r 2 6 Then notice that for weWj and r £ 1Z 

Y D rVi = 



(4.9) 



due to the assumption (4.7). Finally, compute 



E'(Fx;v) = Yl ct>'(Fr)-D r v i = J2<l>'( Fr )- £ ^ = 0. 
(i,i+r) 68 r-e^ i,i+reX 



□ 



In the next subsection we propose a two-dimensional extension of the ECC method (3.3) and 
then show that it is consistent, i.e., that its variation on the uniform deformation y = Fx is zero. 
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4.3 Formulation of ECC 



By analogy with the one-dimensional case, define the continuum contribution of a bond b £ B fully 
contained in fL as 



Cfe(y) := / <KV r y)db (for 6 C O c ) 
Jb 

and hence define the ECC method in 2D: 

E ecc (y) := £ e 6 (y) + £ c b (y). 



(4.10) 



(4.11) 



6eB 



beB 

bcn c 



Note that Q,(y) is well defined since </>(V r y) is piecewise constant along each b C £l c - 



Proposition 4.2. The A/C coupling (4.11) is consistent; i.e., 

(£ ecc )'(Fx; v) = /or a« F G M 2x2 , and u G 

Proof. Since 

c' (M+r) (Fx; u) = <j)'{¥r) ■ 1 V r udb = 0'(Fr) • D r u, = e' (M+r) (Fx; v) 

J (i,i+r) 

the variation of E ecc is zero: 



(4.12) 



(E ccc Y(Fx;v) = 4(fx;v) + £ e^Fsji;) = ^(Far;«) = 0, 
b£n c 

where E"(Fx;i>) = is due to Proposition 



beB 

bcn c 



4.1 



□ 



The expression (4.11) indeed couples the energy of the system in the atomistic and the contin- 
uum region. However, it still remains to be shown that (4.11) can be efficiently computed. The 
following theorem indeed shows this. 



Theorem 4.3. The energy of the A/C coupling defined by (4.11) can be expressed as 



beB 

b<£n c 



where the definition of Cb(y) is extended to b <£_ O c by (4.21) and (4.14) and E c (y) is defined by 



Theorem 4.3 implies that the complexity of computing E ecc (y) scales as the number of bonds 
b G B such that b (£_ O c (which is, up to a constant factor, the same as the number of atoms in fi a ) 
plus the number of triangles in T ■ Remarkably, the formula (4.13) is exactly the same as for the 
one-dimensional model (3.4) (although the respective objects are now in M 2 ). 

The proof of Theorem |4.3| will be given in the next subsection. 
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4.4 Two-dimensional bond-density lemma 



The proof of Theorem 4.3 follows the line of the proof of Proposition |3.3| We first find a two- 



dimensional analogue of the bond-density lemma (Lemma 3.4) and its corollary (3.7). Based on 



them, we extend the definition of Q,(y) to the bonds b <f_ f2 c in a way that (4.13) holds. 

A significant difficulty in extending the one-dimensional results to the two-dimensional case is 
that the bonds are essentially one-dimensional objects; therefore the sum of their characteristic 
functions equals zero almost everywhere in Q c . This means that we have to look for weaker 



analogues of Lemma 3.4 



(4.14) 



To this end, we define the characteristic function for a polygonal set u C M. in the following 
way: 

y (x) = lim ^ nB ^ 
X " [X) £So \B p (x)\ ' 

where B p {x) is the ball with radius p and center x and | • | is the measure of the set. We note 
that (i) the limit w.r.t. p — > in the definition of X ( x ) exists, and (ii) including/excluding the 
boundary of a polygon u (or any part of it) does not change the point values of x ( x )- 
The characteristic function x T ( x ) f° r a triangle T can be visualized as follows: 



X T (x) 



1 
i 

2 
a 
2n 





if x G interior of T, 
if x G edge of T, 

if a; is a vertex of T with angle a, 
otherwise. 



(4.15) 



Note that for the formulation of the method the values of v (x) at the vertices of u will not be 
important. For the characteristic function y (x) thus defined, we have 



X, 



, (x) = y~] X T {x) for all x G R 2 , 



(4.16) 



TeT 



where the identity is strictly pointwise (i.e., not just almost everywhere). 

Using this characteristic function, it is sufficient to obtain the bond-density lemma for a single 
triangle T. 

Lemma 4.4 (bond-density lemma). Let T be a triangle in M 2 whose vertices belong to the lattice 
1? . Then for any r G Z 2 , r ^0, the following identity holds: 



J (i,i+r) 



X T (x)db 



(4.17) 



Proof. Notice that if (4.17) is valid for triangles T\ and T2 (T\ DT2 = 0), then it is valid for T\ UT2. 
Also notice that both sides of (4.17) are invariant w.r.t. translation of T by a vector s G Z 2 and 
reflection around a point p£Z 2 . 

Using these transformations, we can obtain mT (a copy of T stretched by a factor m G Z) as a 
union of m 2 copies of T shifted by different vectors and possibly reflected. This statement can be 
proved by induction: for m = 1 the statement is trivial. Assuming that we have a partition of mT 
by m 2 copies of T, we can complete it to a partition of (m + 1)T by placing m parallelograms (each 
of which consists of two triangles) next to a side of mT and one additional triangle at the corner of 
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Figure 4.2: An illustration of the induction step of a proof that the triangle 3T can be partitioned 
with nine copies of T. We have a partition of 2T with triangles T, T2, T3, T4, and we complete it 
to a partition of 3T by adding two parallelograms (T5 U Tg and T7 U Tg) and one corner triangle Tg. 



that side (see Figure 4.2 for illustration). Thus, we can partition (m+l)T by m 2 + 2m+l = (m+1) 2 
triangles. 

Hence one has 



X T (x)db 



ffl 



l( E / X mT (x)db + 0(m) 

(i,i+r)CmT 

^( E l + 0(m)) 



(i,i+r)CmT 



\ f ImTl + 0(m) ) = |T| + ©(m" 1 ) 



Here the 0{m) terms appear first from neglecting contributions of the bonds near the boundary 
of mT (whose total measure scales as O(m)) and then from estimating the number of points i for 
which (i, i + r) C mT by |mT|, again making the error proportional to the measure of the boundary 
of mT. Letting m — > 00 finally proves (4.17). □ 



Remark 4.1. Unfortunately, the three-dimensional analogue of Lemma 4-4 i- s n °t true. One can 
take a bond direction r = (1,0,0) and consider a tetrahedron T with vertices (1,0,0), (0,1,0), 



(0,0,1), (0,1,1). Then for each bond (i , i + r) the left-hand side of (4.17) will be zero, although the 



right-hand side of (4.17) is 1/6. 



Our next step in deriving an efficient representation of E ecc (y) consists in formulating a weak 
two-dimensional analogue of the identity (3.7). 

Lemma 4.5. Let r G Z 2 , r 7^ 0, y G Vi(7~), and let <f> be a continuous function M 2 — > R. Then the 
following identity holds: 



E / *n <M V ^)db = [ 4>(V r y)dn. 

i+rezaAM+r) 



(4.18) 
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Proof. The proof is based on Lemma 4.4 and identity (4.16). Notice that V r y is constant in each 
triangle. The following formal computation then proves ( 4.18| ): 



E / x nc HV r y)db= E / E^T^ v ^) db 

, ^n,oJ (i,i + r) C ... ^„oJ (i,i+r) rnrT 



' (i,i+r) 



i,i+reZ 2 ' 



>(i,i+r) 



TeT 



E^)| T ( E /. *r db ) 

TeT K Li+reZ^ ^ l+r) J 



E 



TeT 



\T\ 



T 



□ 



As a final step toward deriving an efficient representation of E ecc (y), one can notice that due 
to the assumption (4.8), the sum in the left-hand side of (4.18) can be changed to the sum over 

J. Then 



(i,i + r) £ B, which immediately yields the following corollary. 

Corollary 4.6. Let r £ 1Z, r ^ 0, y £ V±(T), and let f be a continuous function M 2 
the following identity holds: 



E / *n c <M v ^) db = / 4>(V r y)dn. 

■ , _\„ „J (i,i+r) c J 



(i,i+r)eB' 



(4.19) 



4.6 



4.6 



Having established the two-dimensional analogue of partition of unity (Lemma 4.5 and Corollary 
, we can split the continuum energy into individual contributions. For that, apply Corollary 
to the continuum energy (|4.5l): 



E M = E / <Kv r y)dn = E E / x nc <P(v r y)db 

ren^ c rell (i,i+r)eB J ( i > i + r ) 

= E/xn c ^ V ^) db = E c ^)' 



(4.20) 



beB" 



beB 



where we defined 



Cb(y) ■= j^x Qc <P(^ry)db. 
Proof of Theorem \4.3\ Using ( |4.11[ ) and ( |4.20[ ) , compute 

E ecc (y) - E c {y) = E e b (y) + E - E c ^) = E ~ E 



(4.21) 



beB 



tee 



6eB 



fees 

b<£n c 



beB 
b<£Q. c 



which proves (4.13). 



□ 



Remark 4.2 (on an alternative mesh). As a possible modification of the above method, instead 
of requiring that the mesh nodes in the continuum region coincide with the lattice sites, one can 
require that the mesh nodes coincide with lattice half- sites (more precisely, the dual lattice sites), 
i.e., that the nodes of mesh triangles T belong to the lattice Z 2 + (1/2,1/2). One can follow the 
proof of Lemma 4-4 an d show that in this case it is also possible to construct a triangle mT with 
m 2 copies of T reflected around integer points and shifted by integer vectors. 
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4.5 On implementation and relation to QCE 

We can write the continuum energy of ECC using the bond-density lemma (Lemma 4.4) as 

E e c cc (y) ■■= E = E E ^ry\ T )l x T db 

beB TeT beB J(t,*+r) 

bcO, c b(Zfl c 

= E / x T dh) 

TeTren v (i,i+r)eB J ^ l+r ) J 

(i,i+r)cn c (4.22) 

= EE^ v *)(Vi- E / ^ db 

TeTren V (i,i+r)60 Av+r) 

(i,i+r)£Sl c 

= E E ^T,r0(V r y| T ), 



TeTren 



where 



T , r := \T\ - > i x^db 



, a J (i,i+r) 



./ (i,i+r j 

(i,i+r)0f2 c 

is a bond-dependent effective area of T. Notice that £lT,r = \T\ if T is distant enough from the 
interface T. The ECC method can thus be expressed as 

E ecc (y) = E + E E ^T,r<KV r y| T ). (4-23) 

beB TeTren 
be\n c 

The QCE method uses a similar form of the continuum energy: 

E qc (y):=Y,U E e dJ)(y)+ E e (iJ) ( y )) + E^ c E^ v ^)> 

ieZa ^ iex jex ' TeT ren 

j-ien i-jen 

where f2^ c is a particularly defined effective area of T, which also equals |T| if T is distant enough 
from the interface. Thus, the ECC method can be interpreted as a modification of the QCE method 
by allowing the effective areas of the triangles to depend on a bond direction r. 

The implementation of ECC is hence similar to the implementation of the QCE method for two- 
body interaction potentials with the only difference that one needs to compute (or precompute) 
the bond-dependent effective areas £lT,r, which can be done as follows: 

Step 1. Set all &T,r = 1^1- 

Step 2. Loop over all T £ T such that dist(T, T) < R and for each T loop over all bonds b = (i,i+r) E 
B such that b tt c and b intersects with the interior or an edge of T. In the case if b intersects 
with the interior of T, subtract ^ygp from f2;r,r; where |(xi,X2)| is the length of the interval 

between x\ and x%. In the case if b intersects with an edge of T, subtract | from Sl-T.r- 

Note that the factor | is due to X T ( X ) = \ for i G & fl dT. 
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Figure 5.1: An illustration of splitting a bond into several intervals in 2D. The bond b is split 
into three intervals, one with full contribution to the atomistic energy (wi = 1), the other with 
half contribution (w2 = 1/2), and the third one with no contribution (u>3 = 0). The "weights" 



w\,W2,ws are defined in (5.1). 



5 Extension of ACC to 2D 



In section 5.1 we will formulate another two-dimensional coupling method, which we call ACC, 
whose continuum energy is exactly E c {y). The method is a generalization of the respective one- 



dimensional method (3.11). For its two-dimensional generalization we need to define the atomistic 



contribution of a bond coherent with the continuum contribution (4.21). 



Then in section 5.2 we will discuss implementation of ACC. 



5.1 Formulation of ACC 

We start with noticing that Xn (x) is piecewise constant on a given bond b = (i,i + r). In other 
words, there exist = to < t\ < ■ ■ ■ < Im = 1 such that 



!» m := Xn (i + rt )\t^(t t ^ = const for m = 1, . . . ,M. 

lia 1 T t \trri — 1 



(5.1) 



The illustration of splitting a bond b into several intervals is shown in Figure 5.1 Then we define 
the atomistic contribution 



M 



a b(y) = ^2 Wm ( tm ~ 1 



m—l j 



m=l 

and the ACC method 



M 



M 



m—l. 



Yl w m{yi+rt m ~ Vi+rt m -i) I (5.2) 



m=l 



m=l 



£ acc (y):=E a ^) + E c ^)= E^) + £c(y). (5.3) 

beB beB 



beB 

b^n c 



Remark 5.1. If the A/C interface V is convex, then the calculation of a\,(y) can be greatly sim- 
plified, since in that case not more than one section ofb can have a nonzero weight w m , and if the 
weight is not equal to one, then this section of b lies on the interface T. 
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Proposition 5.1. The A/C coupling (5.3) is consistent. 
Proof. Compute the variation of a^{y): 

M i M . 

a' b (Fx; v ) = y~] w m (t m - t m -i) <f/(Fr) ■ ( ^ w m (t m - t m -i) J 

m=l ^ m=l 

M 

= <P\Fr) ■ ^ W m {Vi +r t m - UHrtm-l)- 



-1 A/ 



m=l 



m=l 



Using the fact that V r v (i + ri) = 4j.v(i + rt) (cf. (4.1 )), compute the variation of c b (y): 
c b (Fx;v) = j- X Qc <t>' (Fr) ■ V r vdb 

= J x n (i + rt)(j)'(Vr) ■ V r v(i + rt)dt 

M f tm d 
= Yl / (l-X n (* + rt))0'(Fr)-— v(* + rt))dt 

m =l^*m-l <h 

m=l 

M 

= - w m )4f{?r) • (v i+r t m - ^ +rtm _ 1 ). 

Observe that 

a' b (Fx;t;) + c' b (Fx;v) = 4>{¥r) ■ ^2(v i+r t m - Wi+rt m _i) = 4>'(Fr) • (w i+r - v*) = e' fe (Fx;u). 



m=l 



'■tm ^ 

-r) ■ l ^- v (i + rt))dt 
t i dt 



M 



m=l 



Hence, 



(£ acc )'(Fx; v) = J2 e b (Fx; v) = E'(Fx; v) = 0; 



i.e., the approximation (5.3) is consistent. 



□ 



Remark 5.2. The choice of atomistic contribution (5.2) is not unique. One can, for instance, 
define 



M 



/ Vi+rtm — Ui+rt 

n 



and notice that 



m=l 



M 



a b {Fx;v)= ^ w m (t m - t m -i) $ (Ft 



m=l 



tm tm—1 



m — 1 

tm tm—1 



M 



4>'(Fr) ■ Y w m (v i+rtm - Ui +rtm -i) = a' b (Fx;' 



m=l 
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Figure 5.2: An illustration of three types of bonds that should be considered when implementing 
ACC for convex atomistic region a : (i) bonds that intersect T at one point, (ii) bonds that intersect 
r at two points, and (iii) bonds whose intersection with T is an interval. 



which makes the following alternative version of (|5.3 



E acc (y):=^2~a b (y) + ^2c b (y)= ^ a b (y) + E c (y), (5.4) 



beB beB beB 

b<£n c 



consistent as well. However, in relation to Remark 5.1, one can show that the two methods will 
coincide in the case of a convex atomistic region. 



5.2 Notes on implementation 

In this section we discuss some aspects of implementation of the ACC method (|5.3l) in 2D. 



The ACC energy ( 5.3 ) can be thought of as a sum of continuum energy E c (y) and the respective 
atomistic contributions over the bonds which are not fully inside f2 c (i.e., bonds that are inside fi a 
or crossing the interface T). The continuum energy is treated exactly as in finite elements. We 
assume a triangulation of f2 c whose nodes coincide with lattice sites (this will be referred to as 
nodal atoms) and compute the needed quantities: energy, forces, and stiffness matrix entries; the 
latter is required if Newton-like methods are employed for computing equilibrium of energy. 

The atomistic contributions can be split into two groups, the first one involving only bonds in 
f2 a and the second one involving bonds intersecting with T: 

K cc {y) ■= E a ^ = E <y) + E <y) = : Ki(y) + Ki(y)- ( 5 - 5 ) 

beB beB _ beB 

b<£n c 6cii a &nrv0 

The first sum is nothing but the standard sum of energies of bonds for atoms in S7 a : 



Ki(y)= E e (M)(y)= E <p(\vs-vi\)- 



(tJ)6B (i,j)eB 

(ij)cn a (i,j)cn a 
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The second sum in (5.5) requires extra care. In a general setting, one has to go through all the 
bonds: for each bond find all intersections with the interface T and compute the bond's contribution 
according to formula (5.2). According to Remark 5.1 the implementation can be simplified further 
if the atomistic region Q a is convex (or consists of convex disjoint sets), which is indeed the case 
for most of the simulations of localized defects. In this case all the bonds b such that b D T ^ can 
be further categorized into (i) bonds that intersect T at one point, (ii) bonds that intersect T at 
two points, and (iii) bonds whose intersection with T is an interval (see illustration on Figure 5.2). 

For the type-(i) bonds, we find the point of intersection of a bond b with the interface and 
compute the reconstruction of this point in terms of positions of nodal atoms on the interface. 
Using this reconstruction, we can compute the corresponding contributions to energy, forces, and 
the stiffness matrix. The type-(ii) bonds are treated in the same way, except that we need to 
compute reconstruction of both intersection points. The type- (iii) bonds are treated similarly to 
type (ii) bonds: one should find the reconstruction of two endpoints of an intersecting interval and 
hence compute the contributions to energy, forces, and the stiffness matrix. According to formula 
(5.2), one has a factor of w m = \ for such bonds (since x n | r = 2> (5-1) and ( 4. 15| ) ) . The list of 
coefficients of reconstruction can be precomputed once in the beginning of computation (and once 
per each mesh adaptation if the latter is used). Then, at each iteration, one has to go only through 
this list for assembling the contributions from bonds of type (i), (ii), and (iii). 



6 Numerical experiments 



A number of numerical computations were conducted to illustrate the performance of the proposed 
ECC (cf. ( |4.13[ ) or ( |4.23| )) and ACC (defined in methods in 2D and to compare them with 

the QCE [33] ■ 

In all the numerical examples either the Lennard- Jones potential 



(p(z) 



or the Morse potential 



<p(z) 



-2e 



-2z- 6 + z- 12 



-a(z-l) +e -2a(z-l) 



(6.1) 



(6.2) 



was used. With these potentials, the hexagonal lattice forms a stable equilibrium. 

The test problem was chosen as follows. We took a hexagonal atomic crystal; each side of the 
hexagon contains 129 atoms, and the total number of atoms in the system is 49529. The reference 



atomistic configuration is illustrated in Figure 6.1 Eight atoms has been removed from a perfect 
crystal to form a defect in its center. The atomistic regions formed a smaller hexagon also centered 
at the origin whose side contained K atoms, as illustrated in Figure [l (b) for K = 6. The mesh was 



fully refined near the A/C interface for the ease of implementation of QCE. 

Dirichlet-type boundary conditions were set: we extended the size of the hexagonal region by 
3 atomic spacings and fixed the positions of the three added layers of atoms so that every free 
atom in the system had the full set of neighbors to interact with; this is done in accordance with 

assumptions (4.7) and (4.8). The "external" deformation gradient F = ( M was applied to 



v 0.97 y 

the positions of the added boundary atoms. Such an external compression makes atoms occupy 



the empty lattice sites and form a defect, as illustrated in Figure 6.2 
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Figure 6.1: A reference configuration of the system of 49529 atoms forming a hexagon with the 
side of 129 atoms. The triangulation is shown on the left, and the closeup of the atomistic region 
is shown on the right. Black atoms correspond to degrees of freedom of the system, while white 
atoms are kinematically "slaved" to the black atoms. Eight atoms have been removed to form a 
defect. The mesh is fully refined near the interface. The illustration is for K = 6. 




Figure 6.2: A computed deformation of the atomistic system, closeup of the atomistic region. The 
computation was done with the ECC method for K = 6. 
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In the numerical tests we compute the error of a local minimizer of the coupled A/C methods 
close to a precomputed minimizer of the exact atomistic energy. More precisely, we first compute 
a local minimizer of the exact atomistic energy. We then use it as (i) an initial guess for each 
computation with an A/C method, and (ii) as a reference solution to calculate the error of the A/C 
solution. Of course, in "real" problems we cannot initialize a computation with the exact atomistic 
solution, as it is not known; we do this only to numerically study an approximate solution close to 
the exact one. 

The absolute error of the numerical solution was calculated in the discrete W 1,00 -norm. For the 
test problem used in this work, the W 1,00 -norm of the difference between the reference atomistic 



configuration (Figure 1(b) ) and the exact solution (Figure 6.2) was close to 1; therefore the relative 
error is of the same magnitude as the absolute error. 

A nonlinear conjugate gradient solver with line search |29j was used to find a stable equilibrium 
of an atomistic system. A simple Laplace preconditioner was used to accelerate the convergence. 

The discrete VF 1,00 -norm of the error of an approximate deformation was defined as follows: 
for each triangle formed by the neighboring atoms in the reference configuration we compute the 
Jacobian matrix of the mapping given by the approximate deformations of those atoms, take the 
difference to the corresponding Jacobian matrix of the exact solution, and compute the Frobenius 
norm of this difference. Taking maximum over all triangles of neighboring atoms yields the discrete 
Ty 1 '°°-norm. 



Three tests were conducted: the first one with the Lennard-Jones potential (section 6.1), the 



second one with a slowly decaying Morse potential (section 6.2 ) , and the third one with the Lennard 



Jones potential and a perturbed A/C interface (section 6.3). The results of these computations are 
discussed in section 

6.1 Test with the Lennard— Jones potential 



In the first test case we let atoms interact with the Lennard-Jones potential (6.1) with the cut-off 
distance R = 3.1. We computed the solutions for 5 < K < 50 with the two methods, ECC and 
QCE, and calculated their errors. 



The results of computations are shown in Figure 6.3 We can see that the ECC method does 
converge whereas QCE fails to converge due to ghost forces. The error of QCE shows an initial 
decrease in error for small sizes of the atomistic region (K = 5,6,7) but remains at a level of 
approximately 5 x 10~ 3 as K is further increased. 

The error magnitude of 5 x 1CP 3 is normally acceptable in the engineering applications. Such 
a small error of QCE was due to the second nearest neighbor interaction being relatively weak 
compared to the nearest neighbor interaction for the Lennard-Jones potential. In the next section 
we will see the results for the test case where the second nearest neighbor interaction is considerable. 

6.2 Test with a slowly decaying Morse potential 



In this test case we chose the Morse potential (6.2) with a = 3 and the cut-off distance R = 3.1. 



The strength of such an interaction decays rather slowly, and the ghost force effects are more 

/0.97 \ 

pronounced in this case. The external compression F = I 95 j WaS S6 ^ aS ^ e boundary 



condition. The exact solution for this test is similar to the one shown in Figure 6.2 



CONSISTENT ATOMISTIC/CONTINUUM COUPLING 



26 



error 

0.010 1 i 
0.005 

0.001 
5 x 10" 4 - 

1 x 10" 4 
5 x 10" 5 - 



^oo OOOOO ooooo 
□ 




-% K 



10 15 20 30 50 



Figure 6.3: VF 1,0 °-errors of the computed solutions for the test with the Lennard- Jones potential. 
Computations were done by the two methods, ECC (cf. (4.13) or ( 4.23| ) and QCE [33J, for different 
sizes of the atomistic region (5 < K < 50) and compared with the exact atomistic solution. For a 
small size of the atomistic region (K = 5,6,7), both methods have a comparable error. However, 
for a larger size of the atomistic region (K ^ 8), the ECC method shows a steady convergence 
whereas the error of QCE remains 0(1) due to ghost forces. 
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Figure 6.4: VF 1 ' 00 -errors of the computed solutions for the test with slowly decaying Morse potential. 
Computations were done by the two methods, ECC (cf. (4.13) or ( 4.23| ) and QCE [33J, for different 
sizes of the atomistic region (5 < K < 100) and compared with the exact atomistic solution. The 
ECC method shows steady convergence as the size of the atomistic region increases. In contrast, 
the error of the QCE method is entirely dominated by the ghost force. 
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Figure 6.5: A/C interface of the test with the nonaligned interface (atomistic domain size K = 8). 



The results of computations are shown in Figure 6.4 As can be seen, for such a slowly decaying 
interaction the QCE method exhibits a rather different behavior: its error stayed at the level of 0.05 
no matter how large the atomistic region was. The error of the ECC method, in contrast, steadily 
decayed with increasing K and showed a behavior similar to the test case with the Lennard- Jones 
potential. 



6.3 Test with a nonaligned interface 

In the last test the errors of the ECC and the ACC methods were compared for an aligned and a 



nonaligned A/C interface. The test case is similar to the first one (section 6.1 ), except for a coarser 



mesh in the continuum region (see Figure 6.5). 



The results of computations are shown in Figure 6.6. One can observe that for neither ECC nor 
for ACC the perturbation of the interface affected the error. Also, the error of ECC was slightly 
less for small K as compared to ACC, but for large K the errors were very close. 



7 Discussion 

The main observation that can be made from the results of computations is that the proposed ECC 
and ACC methods exhibit steady convergence as the atomistic domain size is increased and are 
not sensitive to whether the nearest neighbor interaction dominate or whether the A/C interface is 
aligned with the lattice. The solution by the QCE method, in contrast, does not converge due to 
ghost forces. For a rapidly decaying potential, the error of QCE was within a practically acceptable 
limit but was found to be considerably larger for a slowly decaying interaction. 

In the present paper we assumed that both the partition of the material into atomistic and 
continuum regions and the mesh in the continuum region are given. In practice, a good choice of 
the regions and the mesh are often not known a priori, and one has to rely on certain algorithms to 
determine them adaptively (see, for instance, |19U26| ). However, rigorous a posteriori error bounds 
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Figure 6.6: V^ 1,oc -errors of the computed solutions for the test with the Lennard-Jones potential 



and nonaligned interface. Computations were done by the two methods, ECC (cf. (4.13) or (4.23)) 



and ACC (defined in (5.3)), for different sizes of the atomistic region (6 < K < 50) and compared 
with the exact atomistic solution. The squares and circles correspond to the errors of ECC and 
ACC for the aligned interface, while the crosses and pluses correspond to the nonaligned interface. 
Both methods exhibit steady convergence as K is increased, with almost no difference whether the 
interface was aligned or not. 



were proved only in ID for consistent energy-based coupled methods [H [25] . The purely energy- 
based formulation of the proposed method and its consistency may help with deriving rigorous a 
posteriori error bounds in 2D. 

Unfortunately, the bond formulation of the A/C energy developed in this work, according to 
Remark 4.1 cannot be applied directly to couple the exact atomistic model with the Cauchy-Born 
continuum model in three dimensions (3D). Nevertheless, it is possible to formulate an efficient 
numerical algorithm based on the bond formulation of the A/C energy in 3D [27J. 

The major challenge that remains to be addressed for the present method to become competitive 
with existing methods on realistic engineering problems is the formulation of the method for many- 
body potentials, as the two-body potentials are far from covering all existing atomistic models. In 
this regard, we should note that there exists a method of A/C coupling (namely, the geometrically 
consistent scheme of E, Lu, and Yang |8j) which is free from ghost forces under the restriction of 
planar interface with no corners but for general potentials. 

It would also be interesting to extend the present method to coupling an atomistic model with a 
nondiscretized continuum model in 2D and 3D. In that case one can discretize the continuum model 
with higher-degree polynomials and obtain an increased accuracy of the overall method. Another 
extension which would be useful for applications is to formulate the proposed method for complex 
lattices (for a possible application, see [H] for two-dimensional models of metallic alloys described 
in terms of two-body potentials). 



8 Conclusion 



We considered the problem of consistent coupling of atomistic and continuum models of materials, 
limited to the case of a two-body potential and one or two spatial dimensions. We proposed two 
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versions of energy-based coupling which are consistent (i.e., do not suffer from ghost forces). The 
coupling is based on judiciously defining the contributions of the atomistic bonds to the discrete 
and the continuum potential energies. The same coupling in ID has been independently proposed 
and analyzed in [IB]. The latest works conducted on the proposed coupling include its numerical 
analysis in 2D [24J and extension to 3D [27J. 
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